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A new algorithm has been developed for the estimation of tropospheric microwave 
path delays from water vapor radiometer ( WVR } data, which does not require site and 
weather dependent empirical parameters to produce high accuracy . Instead of taking the 
conventional linear approach, the new algorithm first uses the observables with an emis- 
sion model to determine an approximate form of the vertical water vapor distribution 
which is then explicitly integrated to estimate wet path delays, in a second step . The 
intrinsic accuracy of this algorithm has been examined for two channel WVR data using 
path delays and simulated observables computed from archived radiosonde data . It is 
found that annual RMS errors for a wide range of sites are in the range from 1.3 mm to 
2.3 mm, in the absence of clouds. This is comparable to the best overall accuracy obtain- 
able from conventional linear algorithms, which must be tailored to site and weather 
conditions using large radiosonde data bases. The accuracy and flexibility of the new 
algorithm are indications that it may be a good candidate for almost all WVR data 
interpretation. 


I. Introduction 

Very Long Baseline Interferometry (VLBI) and Global 
Positioning System (GPS) applications, as well as other fields, 
require calibration of microwave delays caused by water 
vapor in the troposphere. Two and three channel water vapor 
radiometers (WVRs) are increasingly seen as the best solution 
to this problem. Unfortunately, WVRs do not directly measure 
delay; they actually measure antenna temperatures in fre- 
quency channels close to the 22.2 GHz resonance of water. 
The WVR observables must then be related to the desired wet 
path delays by a data interpretation algorithm. In order to 


most effectively use the WVR data, one must, therefore, 
understand and quantify the level of error in the delay estima- 
tion contributed by the algorithm itself. 

The most popular algorithm formulation, resulting from the 
conventional approach to delay estimation, assumes that the 
microwave delay may be written as a linear expansion in the 
total optical depths of the troposphere at the WVR channel 
frequencies (e.g., Refs. 1, 3, and 6). The optical depth at any 
point s , for a frequency is defined by the integral along the 
line-of-sight in the atmosphere, 
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The total optical depth to be used in the delay expansion is, 
therefore, just r(°°, Vj). The symbol a represents the extinction 
coefficient, also called opacity (e.g., Ref. 5, page 23), a com- 
plicated function of temperature, pressure, frequency, and the 
densities of water vapor and liquid, that is defined by an atmo- 
spheric emission model (e.g.. Refs. 9 and 10). In the linear 
algorithms, the r(°°, v ( ) are usually estimated from the antenna 
temperatures after assuming effective radiating temperatures 
for the tropospheric water vapor. The coefficients in the 
linear delay expansion, sometimes called retrieval coefficients 
(see Ref. 3), are determined empirically, usually using a set ot 
wet delays and corresponding simulated WVR antenna tem- 
peratures computed from archived radiosonde data. Once the 
delay coefficients have been determined by regression analysis 
of the simulated observables, the WVR user has a very simple, 
easily applied, expression for wet delay. 


There are also undesirable features of linear delay algo- 
rithms, as with most other algorithms. In particular, the linear 
coefficients must be empirically tailored to each observing 
site and set of observing conditions in a process sometimes 
called stratification, in order to achieve high delay estimation 
accuracy (Ref. 3). The need for this arises from the fact that 
the WVR antenna temperatures and surface meteorology data 
do not uniquely define the implicit parameters of the delay 
estimation problem. Specifically, the WVR antenna tempera- 
tures, T a (y^ are imperfect measurements of the radiation 
brightness temperatures, which are given by 


w= 





dr T k (r) exp(-r) + T bg exp(-T(°°, vj) 


or, equivalently, 


T b (Vi) = I ds a (s, v t ) T k (s ) exp(-r(s, uj)+ T bg exp v,)) 

Jo 

where r is defined above, T k is the kinetic temperature, and 
T bg is the background temperature. Because of the depen- 
dence on r and a, the WVR observables also depend on the 
distributions of temperature, pressure, and water liquid, as 
well as the actual amount of water vapor along the line-of- 
sight. The conventionally defined wet path delay has a simpler 
form, 


R 


wet 



ds p v (s)/T k (s) 


but also depends on the distribution of water vapor density, 
p (s), and the temperature profile (e.g.. Ref. 2). One can see 
from these expressions that only a knowledge of the line-of- 
sight temperature, pressure, water vapor, and water liquid 
distributions can uniquely specify both brightness tempera- 
tures and wet path delay. Since this information is not directly 
available in observable quantities, the past values of the required 
data under specific conditions are used to improve the defini- 
tions of the linear delay coefficients. 

Stratification of the coefficients requires that a large radio- 
sonde data base be available for the empirical estimation of 
each coefficient set. The first difficulty that this presents is 
that of simply matching coefficients to individual sites. Since 
radiosondes are flown at relatively few locations which are 
usually near population centers and away from most radio 
observatories, there is nearly always some uncertainty in the 
association of a data base with any specific site of interest. 
In addition, even with the currently discussed computation of 
delay coefficients for specific seasons, weather conditions, 
and sites, use of the linear algorithm may become awkward 
when processing typical VLBI and GPS experiments. Finally, 
the delay estimation accuracy of the linear algorithms cannot 
be improved indefinitely by simply increasing the level of 
stratification, because of the finite number of meaningful 
observable criteria for discrimination. 

In order to avoid the disadvantages of linear algorithms, a 
new algorithm, the profile algorithm, has been developed from 
an entirely different approach. The new algorithm’s formula- 
tion of the delay estimation problem relies on the direct 
application of the known physics of the problem to make 
better use of observables, rather than on large data bases to 
adjust linear delay coefficients for specific conditions. The 
result is more complex than the linear algorithms, but has 
considerable advantages. Most notably, it offers complete 
independence from site and condition-specific empirical 
parameters, and their inherent problems, while producing 
delay estimation accuracy that is as good as that of the better 
stratified linear algorithms. 


II. Description of the Profile Algorithm 

The basic premise of the profile algorithm is that an actual 
solution to a simplified version of the nonlinear delay problem 
can be generated numerically from each WVR observation and 
the corresponding surface meteorology. This leads to a formu- 
lation and an execution very unlike those for the linear algo- 
rithms. In fact, the profile algorithm makes no attempt to 
impose a direct relationship, linear or otherwise, between 
WVR antenna temperatures and wet path delays. Instead, 
simple vertical distributions of water vapor and liquid are 
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adjusted until they produce, through an atmospheric emission 
model (e.g., Ref. 10), brightness temperatures in agreement 
with the WVR antenna temperatures. The delay estimates 
themselves are produced in a second step by integrating along 
the line-of-sight using the fitted water vapor distributions. 

It is assumed that the vertical distribution of the relative 
humidity can be adequately represented by a two piece linear 
function. At surface level, the relative humidity is naturally 
given by surface meteorology, and at altitude of 10 km, so 
little water vapor can exist that it can safely be taken as 
zero. The relative humidity profile is then approximated by 
connecting the two end points linearly to the relative humidity 
at an intermediate altitude of 3 km. It is the value of the 
relative humidity at this point that the algorithm freely 
adjusts in order to completely specify the vapor profile. An 
example of such a simplified water vapor profile used by the 
profile algorithm is illustrated on the left side of Fig. 1. The 
corresponding water vapor density profile, defined by the 
sample relative humidity distribution, is shown in the right 
portion of the figure. 

The vertical pressure and temperature profiles must also be 
known, or assumed, in order to specify the water vapor 
density and opacity from the algorithm’s relative humidity 
profile. The pressure can be well represented by an exponen- 
tial decay from the surface value. The temperature profile for 
the profile algorithm is produced by first creating a nominal 
temperature profile from data in Tables 5.1 and 5.2 of the 
Standard Atmosphere (Ref. 12), interpolated to the latitude 
and altitude of the site, and to the time of year of the observa- 
tion. The surface temperature implied by the nominal profile 
is then subtracted from the measured surface temperature 
and multiplied by exp (-h/H), where h is the altitude, and H is 
a constant scale height, to produce an altitude-dependent 
temperature correction. The temperature profile used by the 
algorithm is then formed by adding the nominal profile and 
the correction. It, therefore, matches the measured surface 
temperature, and exponentially approaches the nominal 
temperature distribution as the altitude increases. An examina- 
tion of temperature profiles from radiosondes at a number of 
locations suggests that the best choice of scale height for the 
decay of the surface temperature adjustment is about 2 km. 
An example of such an adjusted profile is shown in Fig. 2. 
Note that the site altitude and latitude, needed to produce the 
nominal temperature profiles, are the only site-specific param- 
eters required by the profile algorithm. 

In reality, the actual vertical water vapor distributions 
deviate from the assumed form of the profiles; however, the 
impact of these differences is quite small, as the next section 
of this article demonstrates. One reason for this is that the 
channel frequencies of the WVRs are usually selected so as to 


minimize the influence of changes in the vertical water vapor 
profiles on estimated delay (Refs. 7 and 1 1). Another is that in 
integrating the vapor distribution to find the wet delay, high 
frequency spatial variations tend to average out, making a 
smooth function an adequate representation. Still, dependence 
on the vertical structure cannot be completely eliminated, so 
the fact that the new algorithm provides even an approximate 
adjustment to the effective height of the water vapor gives 
some improvement over the linear algorithms. 

In the simulations that are described in the next section, it 
is assumed that the liquid opacity is negligible, but in actual 
use the delay estimation algorithm must accurately account 
for the continuum emission from water droplets in clouds in 
order to produce high quality delay estimates. The liquid 
opacity is not a function of pressure, and only weakly depen- 
dent on temperature (e.g., Ref. 9), which means that the 
liquid contribution to the WVR antenna temperatures is 
relatively insensitive to the actual form of the vertical distri- 
bution of droplets. Therefore, the simple assumption that the 
water liquid density is proportional to the saturation water 
vapor density at any point has been made, and the algorithm 
is allowed to fit for the constant of proportionality. The effect 
of this assumption will be examined in work, unfinished at 
the time of this writing, dealing with algorithm performance 
in the presence of clouds. 

III. Accuracy of the Profile Algorithm 

All algorithms for the estimation of wet path delays from 
WVR data contribute some error to the final estimates. This 
error is intrinsic to the algorithm, and as such can be separated 
from other error sources. The best way to isolate and quantify 
the algorithm error is by generating simulated WVR antenna 
temperatures and corresponding wet path delays from a set 
of vertical profiles of water density, temperature, and pres- 
sure. This is a process similar to that used to generate empirical 
coefficients for linear algorithms, as described by Gary et al. 
(Ref. 3), Resch (Ref. 6), and others. Such a simulation has 
been done for the new profile algorithm, and the results are 
described in this section. 

Radiosonde data are taken on a regular basis at a large 
number of sites over the world, and can be used as the basis 
of the required simulation. These data typically consist of 
temperature and relative humidity on a vertical grid of pres- 
sure, which can be converted to altitude. For our simulation, 
data were obtained spanning a full year for 16 sites of interest. 
Some sites were selected from the available radiosonde archive 
because they were closest to prominent VLBI observatories, 
others were chosen to overlap with other delay retrieval 
accuracy studies, and the rest were picked to produce a good 
cross-section of locations in the continental United States. 
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Although 13 of the data sets contained two radiosonde pro- 
files for most days, Edwards AFB, Madrid, and Wagga-Wagga 
had less than half that number available. These were, nonethe- 
less, included because of their proximity to Goldstone, DSS 63, 
and DSS 43, respectively. Radiosonde profiles which had 
missing or blatantly erroneous data below 5 km, or which 
terminated at anomalously low altitudes, were removed 
by software to avoid corruption of the results. This left an 
average of 645 profiles for all but the three previously men- 
tioned sites, which averaged 260. Examination of the re- 
jected data revealed no selection effects that could bias the 
simulation. 

Next, each profile was used to compute a zenith delay, 
due to water vapor, and the corresponding antenna tempera- 
tures that would be observed by an ideal noise-free WVR with 
frequency channels at 20.7 GHz and 31.4 GHz. This par- 
ticular WVR configuration was selected because it corresponds 
to the majority of existing water vapor radiometers. The 
emission model used to calculate the WVR observables is 
essentially that of Waters (Ref. 10) with oxygen opacities 
given by Rozenkranz (Ref. 8). Because the same emission 
model was used for both the simulation and the delay re- 
trieval, all emission model errors should cancel in the simula- 
tion. Of course, these errors will contribute to the ultimate 
delay estimation accuracy, but they are the same for all algo- 
rithms. No cloud contribution to the continuum emission was 
used, at this time, because no information concerning liquid 
content is present in the radiosonde profiles, and a suitable 
cloud model had not yet been implemented. 

Finally, the simulated WVR antenna temperatures and sur- 
face meteorology data were used as input to the profile 
algorithm to produce the delay estimates. The level of algo- 
rithm error is then given by the difference between the retrieved 
delay estimates and the corresponding delays computed 
directly from the radiosonde profiles. These results for all 
16 sites appear in Table 1. The last column displays the root- 
mean-square algorithm error for the entire year at each loca- 
tion. These range from 1.3 mm for Wagga-Wagga up to 2.3 mm 
for Munich, and have an average value of 1 .8 mm. The annual 
average zenith delay and the maximum error in estimated 
delay are also displayed in the table. 

Resch presented several variations on the linear delay 
retrieval algorithm in his 1983 paper (Ref. 6). The most 
popular, and most accurate, makes use of surface meteorology, 
and fits the standard form of the linear algorithm given below. 

^wet = c O + a 0 + a \ T (°°- "l) + a 2 T (°°- V t) 

R wet is the wet path delay estimate, the a i and c Q are empiri- 
cally determined constants, and jj. is the cosine of the zenith 


angle of the observation. The r(°°, v t ) are optical depth esti- 
mates, corresponding to the WVR frequency channels, com- 
puted from 

T(°°, V t ) = -In [(T eff (V,.) - T a (y$)KT eff (yd - 2.9)] 

The r eff (^) in the above equation are empirically determined 
effective radiating temperatures, computed by fitting to 
simulated WVR observations, and the T a (y f .) are the measured 
antenna temperatures. 

For this algorithm, Resch determined best-fits to site- 
specific linear delay coefficients using 90 radiosonde profiles 
to cover a one-year period for each of five different locations. 
He then presented the root-mean-square errors of those fits 
to the data bases under the same assumptions of no liquid 
opacity and no WVR noise that were used in the simulations 
here. The results in this article may, therefore, be compared 
directly to his for these five sites, as has been done in Table 2. 
Although the errors for the individual sites are generally 
different, the average RMS error, at 2 mm, is not significantly 
different for the stratified linear algorithm and the profile 
algorithm with no empirically determined site-specific param- 
eters, at all. 

IV. Discussion and Conclusions 

Because of its need to perform emission model calculations, 
the profile algorithm is required to do a relatively large amount 
of computation in order to estimate delays from WVR obser- 
vations. Linear algorithms, on the other hand, require exten- 
sive calculations to initially produce linear delay coefficients, 
but thereafter need very little computation to produce delay 
estimates. This is not a serious drawback for the new algorithm, 
however, since even the microcomputers that run the more 
recent WVRs, such as Janssen’s J01 (Ref. 4), are fully capable 
of running real-time software implementing the profile algo- 
rithm. The emission model calculations also produce a definite 
advantage in that the tropospheric emission characteristics 
are imperfectly known, and subject to revision. Any emission 
model changes made to the profile algorithm code can imme- 
diately be used for reducing data without delay. In contrast, 
linear algorithms require a repetition of the lengthy simula- 
tions, reprocessing all of the radiosonde data originally used to 
generate the linear delay coefficients, before new delay esti- 
mates may be produced. 

There may well be instances when WVR users would prefer 
to spend only the smaller amount of computation time needed 
by the linear algorithms, without accepting their less desirable 
features. For such situations, suitable hybrid schemes can be 
devised. For example, the profile algorithm is particularly 
well-suited for use with tipping-curve reduction programs; a 
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very simple extension would allow it to use the tipping-curve 
data to simultaneously produce linear delay coefficients, 
which could then be used for delay estimation until the next 
invocation of the profile algorithm. There are, of course, 
many other possible ways of making the best use of both types 
of algorithms in combination. 

The work presented here demonstrates the validity of the 
profile algorithm and, since many other related formulations 
are possible, the approach to the WVR delay estimation 
problem that the new algorithm represents. The most serious 
question still to be addressed is the profile algorithm accuracy 


in the presence of clouds. Preliminary indications from work 
now in progress are that, while cloudy weather does decrease 
the delay estimation accuracy, the overall algorithm accuracy 
is not degraded by more than 15 percent by the addition of 
liquid opacity to the simulations. Even without the liquid 
opacity information, the flexibility of the new algorithm, and 
its apparent ability to use the same procedures and parameters 
to reduce WVR data from any site under a wide range of 
weather conditions, while still obtaining high delay estimation 
accuracy, may make it a good candidate for use with almost all 
WVR applications. Additional work to thoroughly investigate 
its properties is certainly justified. 
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Table 1. Profile algorithm errors for 1979 


Site 

Average delay, cm 

Maximum error, cm 

RMS error, cm 

Albuquerque 

7.4 

0.55 

0.17 

Apalachicola 

18.0 

0.85 

0.19 

Boise 

7.4 

0.52 

0.16 

Dayton 

12.3 

0.90 

0.20 

Denver 

6.6 

0.62 

0.14 

Edwards AFB 

6.6 

0.51 

0.17 

El Paso 

8.7 

0.92 

0.17 

Goteborg 

8.0 

0.77 

0.16 

Madrid 

10.4 

0.95 

0.19 

Munchen 

10.0 

0.81 

0.23 

Oakland 

9.0 

0.64 

0.20 

Oklahoma City 

12.5 

0.93 

0.20 

Pittsburgh 

11.0 

0.77 

0.18 

Portland (ME) 

10.3 

0.97 

0.22 

San Diego 

9.0 

0.84 

0.22 

Wagga-Wagga 

9.1 

0.41 

0.13 


Minimum RMS error 

0.13 


Maximum RMS error 

0.23 


Average RMS error 

0.18 


Table 2. Comparison of errors from profile and Resch algorithms 


Site Average delay, cm 

Profile RMS error, cm 

Resch RMS error, cm 

El Paso 8.7 

0.17 

0.14 

Oakland 9.0 

0.20 

0.27 

Pittsburgh 11.0 

0.18 

0.16 

Portland (ME) 10.3 

0.22 

0.18 

San Diego 9.0 

0.22 

0.30 

Minimum RMS error 

0.17 

0.14 

Maximum RMS error 

0.22 

0.30 

Average RMS error 

0.20 

0.21 


RELATIVE HUMIDITY 


VAPOR DENSITY, g/cm° X 10"' 


Fig. 1 . An example of the assumed relative humidity profile appears on the left. The dashed line at 3 km represents 
the value of the fitting parameter. The water vapor density profile corresponding to the relative humidity is shown on 
the right. 
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Fig. 2. An example of a nominal temperature profile derived from 
data in the Standard Atmosphere is represented by the dashed line. 
The solid line shows the temperature profile after it has been forced 
to match the measured surface temperature. 
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